home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / gfsr4.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  5.8 KB  |  173 lines
  1. /* This program is free software; you can redistribute it and/or
  2.    modify it under the terms of the GNU General Public License as
  3.    published by the Free Software Foundation; either version 2 of the
  4.    License, or (at your option) any later version.
  5.  
  6.    This program is distributed in the hope that it will be useful, but
  7.    WITHOUT ANY WARRANTY; without even the implied warranty of
  8.    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  9.    General Public License for more details.  You should have received
  10.    a copy of the GNU General Public License along with this program;
  11.    if not, write to the Free Foundation, Inc., 59 Temple Place, Suite
  12.    330, Boston, MA 02111-1307 USA
  13.  
  14.    From Robert M. Ziff, "Four-tap shift-register-sequence
  15.    random-number generators," Computers in Physics 12(4), Jul/Aug
  16.    1998, pp 385-392.  A generalized feedback shift-register (GFSR)
  17.    is basically an xor-sum of particular past lagged values.  A 
  18.    four-tap register looks like:
  19.       ra[nd] = ra[nd-A] ^ ra[nd-B] ^ ra[nd-C] ^ ra[nd-D]
  20.    
  21.    Ziff notes that "it is now widely known" that two-tap registers
  22.    have serious flaws, the most obvious one being the three-point
  23.    correlation that comes from the defn of the generator.  Nice
  24.    mathematical properties can be derived for GFSR's, and numerics
  25.    bears out the claim that 4-tap GFSR's with appropriately chosen
  26.    offsets are as random as can be measured, using the author's test.
  27.  
  28.    If the offsets are appropriately chosen (such the one ones in
  29.    this implementatin), then the sequence is said to be maximal.
  30.    I'm not sure what that means, but I would guess that means all
  31.    states are part of the same cycle, which would mean that the
  32.    period for this generator is astronomical; it is
  33.      (2^K)^D ~ 10^93334
  34.    where K=32 is the number of bits in the word, and D is the longest
  35.    lag.  This would also mean that any one random number could 
  36.    easily be zero; ie 0 <= ra[] < 2^32.
  37.  
  38.    Ziff doesn't say so, but it seems to me that the bits are
  39.    completely independent here, so one could use this as an efficient
  40.    bit generator; each number supplying 32 random bits.
  41.  
  42.    This implementation uses the values suggested the the author's
  43.    example on p392, but altered to fit the GSL framework.  The "state"
  44.    is 2^14 longs, or 64Kbytes; 2^14 is the smallest power of two that
  45.    is larger than D, the largest offset.  We really only need a state
  46.    with the last D values, but by going to a power of two, we can do a
  47.    masking operation instead of a modulo, and this is presumably
  48.    faster, though I haven't actually tried it.  The article actually
  49.    suggested a short/fast hack:
  50.  
  51.    #define RandomInteger (++nd, ra[nd&M]=ra[(nd-A)&M]\
  52.                           ^ra[(nd-B)&M]^ra[(nd-C)&M]^ra[(nd-D)&M])
  53.  
  54.    so that (as long as you've defined nd,ra[M+1]), then you ca do things
  55.    like: 'if (RandomInteger < p) {...}'.
  56.  
  57.    Note that n&M varies from 0 to M, *including* M, so that the
  58.    array has to be of size M+1.  Since M+1 is a power of two, n&M
  59.    is a potentially quicker implementation of the equivalent n%(M+1).
  60.  
  61.    This implementation copyright (C) 1998 James Theiler, based on
  62.    the example mt.c in the GSL, as implemented by Brian Gough.
  63. */
  64.  
  65. #include <config.h>
  66. #include <stdlib.h>
  67. #include <gsl/gsl_rng.h>
  68.  
  69. static inline unsigned long int gfsr4_get (void *vstate);
  70. static double gfsr4_get_double (void *vstate);
  71. static void gfsr4_set (void *state, unsigned long int s);
  72.  
  73. /* Magic numbers */
  74. #define A 471
  75. #define B 1586
  76. #define C 6988
  77. #define D 9689
  78. #define M 16383 /* = 2^14-1 */
  79. /* #define M 0x0003fff */
  80.  
  81. typedef struct
  82.   {
  83.     int nd;
  84.     unsigned long ra[M+1];
  85.   }
  86. gfsr4_state_t;
  87.  
  88. static inline unsigned long
  89. gfsr4_get (void *vstate)
  90. {
  91.   gfsr4_state_t *state = (gfsr4_state_t *) vstate;
  92.  
  93.   state->nd = ((state->nd)+1)&M;
  94.   return state->ra[(state->nd)] =
  95.       state->ra[((state->nd)+(M+1-A))&M]^
  96.       state->ra[((state->nd)+(M+1-B))&M]^
  97.       state->ra[((state->nd)+(M+1-C))&M]^
  98.       state->ra[((state->nd)+(M+1-D))&M];
  99.   
  100. }
  101.  
  102. static double
  103. gfsr4_get_double (void * vstate)
  104. {
  105.   return gfsr4_get (vstate) / 4294967296.0 ;
  106. }
  107.  
  108. static void
  109. gfsr4_set (void *vstate, unsigned long int s)
  110. {
  111.   gfsr4_state_t *state = (gfsr4_state_t *) vstate;
  112.   int i, j;
  113.   /* Masks for turning on the diagonal bit and turning off the
  114.      leftmost bits */
  115.   unsigned long int msb = 0x80000000UL;
  116.   unsigned long int mask = 0xffffffffUL;
  117.  
  118.   if (s == 0)
  119.     s = 4357;    /* the default seed is 4357 */
  120.  
  121.   /* We use the congruence s_{n+1} = (69069*s_n) mod 2^32 to
  122.      initialize the state. This works because ANSI-C unsigned long
  123.      integer arithmetic is automatically modulo 2^32 (or a higher
  124.      power of two), so we can safely ignore overflow. */
  125.  
  126. #define LCG(n) ((69069 * n) & 0xffffffffUL)
  127.  
  128.   /* Brian Gough suggests this to avoid low-order bit correlations */
  129.   for (i = 0; i <= M; i++)
  130.     {
  131.       unsigned long t = 0 ;
  132.       unsigned long bit = msb ;
  133.       for (j = 0; j < 32; j++)
  134.     {
  135.       s = LCG(s) ;
  136.       if (s & msb) 
  137.         t |= bit ;
  138.       bit >>= 1 ;
  139.     }
  140.       state->ra[i] = t ;
  141.     }
  142.  
  143.   /* Perform the "orthogonalization" of the matrix */
  144.   /* Based on the orthogonalization used in r250, as suggested initially
  145.    * by Kirkpatrick and Stoll, and pointed out to me by Brian Gough
  146.    */
  147.   for (i=0; i<32; ++i) {
  148.       int k=7+i*3;
  149.       state->ra[k] &= mask;    /* Turn off bits left of the diagonal */
  150.       state->ra[k] |= msb;    /* Turn on the diagonal bit           */
  151.       mask >>= 1;
  152.       msb >>= 1;
  153.   }
  154.  
  155.   state->nd = i;
  156. }
  157.  
  158. static const gsl_rng_type gfsr4_type =
  159. {"gfsr4",            /* name */
  160.  0xffffffffUL,            /* RAND_MAX  */
  161.  0,                    /* RAND_MIN  */
  162.  sizeof (gfsr4_state_t),
  163.  &gfsr4_set,
  164.  &gfsr4_get,
  165.  &gfsr4_get_double};
  166.  
  167. const gsl_rng_type *gsl_rng_gfsr4 = &gfsr4_type;
  168.  
  169.  
  170.  
  171.  
  172.  
  173.